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HIGHLIGHTS 


• Melting of a nuclear fuel rod under the accident conditions is analyzed. 

• Enthalpy formulation with explicit finite difference discretization is adopted. 

• Temperature histories and phase change interface positions are determined. 

• The heat transfer coefficient and heat generation rate are two predominant factors. 

• Classical lumped model fails to give accurate results for high heat transfer coefficient. 
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Melting of nuclear fuel rods is a crucial issue that must be addressed when conducting simulation analyses 
of severe accidents in nuclear reactors. The present work aims to develop a mathematical model and carry 
out numerical simulation of one-dimensional heat conduction with phase change in a nuclear fuel rod 
including two regions: uranium dioxide fuel pellet and zirconium alloy cladding. The temperature- 
dependent specific heats, thermal conductivities and densities are considered in the present investiga- 
tion. The enthalpy formulation of the melting process in the nuclear rod during full reactor power oper- 
ation is solved by using a finite-difference method. The thermal analysis is divided into four different 
phases: (1) transient phase before the cladding melting, (2) cladding melting phase, (3) transient phase 
before the fuel melting and (4) fuel melting phase. Then, the effect of heat transfer coefficient between 
coolant and fuel rod on the temperature histories of the fuel pellet and the cladding is investigated. Finally, 
the melting process of a nuclear rod with decay reactor power after shutdown is also simulated. 

© 2014 Elsevier Ltd. All rights reserved. 


1. Introduction 

Core meltdown, which may happen due to the inadequate 
cooling of the nuclear fuel rods, is a serious accident scenario in 
nuclear reactors, such as the Three Mile Island, the Chernobyl and 
the Fukushima nuclear accidents ]JJ. A core melt accident occurs 
when the heat generated by the nuclear reactor exceeds the heat 
removed by the reactor cooling systems to the point where at least 
one nuclear fuel element exceeds its melting point. In a pressurized 
water reactor (PWR) or boiling water reactor (BWR), the melting of 
fuel rods may release radioactive material into the cooling water 
and cause the zirconium cladding to react with water to generate 
explosive hydrogen. 

Many analytical and numerical approaches have been proposed 
to obtain the transient heat conduction response of the nuclear fuel 
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rod. Based on the assumption that the heat transfer resistance 
between the fuel and clad remains constant for both steady-state 
and transient periods, Chen et al. [2] developed a simple conduc- 
tion model with phase change for the transient analysis of a reactor 
fuel rod based on average properties and lumped parameter tech- 
niques, where the analysis is subdivided into four phases: transient 
phase before the clad melting, clad melting phase, transient phase 
before the fuel melting and fuel melting phase. Ghiaasiaan et al. [3] 
obtained the solution of the governing equations for the radial heat 
conduction in the fuel pellet and cladding by the lumped capaci- 
tance method, the method of lines and the integral method, 
respectively. Regis et al. £41 analyzed the transient heat conduction 
in a nuclear fuel rod by an improved lumped parameter approach, 
where Hi o and Ho,o Hermite approximations were applied for the 
averaged temperatures and heat fluxes. Su and Cotta [5] presented 
a higher order lumped parameter formulation for simplified LWR 
(light water reactor) thermohydraulic analysis that the two-sided 
corrected trapezoidal rule ( H ^ approximation) was used in the 
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Nomenclature 

c p specific heat, J/(kg • I<) 

h convective heat transfer coefficient, W/(m 2 • I<) 

H enthalpy, J/kg 

k thermal conductivity, W/(m • I<) 

L latent heat of melting, J/kg 

g volumetric heat generation, W/m 3 

r space coordinate, m 

R radius, m 

s solid-liquid interface location 

t time, s 

T temperature, I< 

Greek letters 
p density, kg/m 3 

Subscripts 
c cladding 

f fuel 

i node 

1 liquid phase 

m melting phase 

s solid phase 


averaged temperature integrals for the fuel and the cladding, and 
plain trapezoidal rule (H 0) o approximation) was used in the aver- 
aged heat flux integrals. Pontedeiro et al. [6] proposed an improved 
lumped-differential formulation for one-dimensional transient 
heat conduction in a heat-generating cylinder with temperature- 
dependent thermophysical properties typical of high burn-up nu- 
clear fuel rods. To eliminate the paradox of an infinite thermal wave 
speed, Espinosa-Paredes and Espinosa-Martinez [7] examined the 
applicability of a fuel rod mathematical model based on Non- 
Fourier transient heat conduction as constitutive law for the Light 
Water Reactors transient analysis (LWRs). Kudryashov et al. [8] 
studied the temperature distribution in the nuclear fuel rod of 
the reactor taking the influence of the high burn-up into account, 
and found the solution for the temperature distribution in the 
nuclear fuel rod in the analytical form for the stationary behaviour 
of the reactor. To overcome the limitation of the separation of 
variables method, Singh et al. [9J presented the finite integral 
transform method to solve the asymmetric heat conduction prob- 
lem in a multilayer annulus with time-dependent boundary con- 
ditions and/or heat sources. 

The melting behaviour of fuel rods during an accident is basi- 
cally determined by solving the heat conduction equation 10]. The 
typical approaches to investigate the phase change problems 
include finite-difference (FD) scheme [111 , heat balance integral 
method [12], moving mesh approach Q3J, variable iteration 
method [141 , lumped parameter model [151 , etc. For example, 
Voller and Cross [111 presented an explicit FD scheme for the 
enthalpy formulation of one-dimensional freezing problems, and 
also developed an implicit FD algorithm of the Stefan problems 
with internal heat generation and melting temperature defined by 
a mushy region, simultaneously. Cheung et al. [161 studied 
numerically the process of freezing and melting occurring in a heat- 
generating slab bounded by two semi-infinite cold walls, and 
solved the dimensionless governing equations using the method of 
collocation with Hermite splines as approximating functions and 
Gaussian quadrature points as the collocation points. Using a quasi- 
steady approximation, Jiji and Gaye [171 examined analytically one- 
dimensional solidification and melting of a slab with uniform 


volumetric energy generation, and applied the results to two ex- 
amples: solidification of a nuclear material and melting of ice. With 
the perturbation technique, Yu et al. [181 studied the planar solid- 
ification with time-dependent heat generation in a semi-infinite 
plane, where the results of the proposed method were validated 
by the exact solution of the classical Stefan problem without heat 
generation. To assess the recriticality potential for a liquid metal 
fast breeder reactor (LMFBR) following a core disruptive accident, 
El-Genk and Cronenberg [191 studied the freezing phenomena of 
molten fuel on a cold structure, and used the successive approxi- 
mation technique to obtain a solution to the non-linear freezing 
problem. The effects of heat generation, viscous heat dissipation, 
temperature-dependent thermophysical properties and a convec- 
tive boundary condition at the solidification front have been 
incorporated into the an analytical formulation. Crepeau et al. [201 
derived approximation governing equations of the interface loca- 
tion in one-dimensional PCM structure with internally generated 
heat in cylindrical, spherical, plane wall and semi-infinite geome- 
tries. Extending previous work, Crepeau et al. [211 investigated the 
solid-liquid phase change driven by volumetric energy generation 
in a vertical cylinder by a quasi-static, approximate analytical so- 
lution, and studied numerically the effect of convection within the 
liquid region. Crepeau and Siahpush [ 221 presented a quasi-static 
analytical solution of melting process in a cylinder with volu- 
metric heat generation, which was valid for the Stefan number less 
than one. 

This work aims at analyzing heat conduction for the accident 
scenarios in the core of a pressurized water reactor (PWR), specially 
the thermal behaviours of uranium dioxide pellet and zirconium 
alloy tubing under high temperature conditions that can cause the 
core to melt. Based on the enthalpy formulation, the governing 
equation is discretized by explicit FD methods with the forward and 
central differential schemes for time and space derivatives, 
respectively. The transient heat conduction of the fuel pellet with 
internal heat generation is considered, while there is no heat 
generation in the cladding. The temperature-varying thermophys- 
ical properties such as thermal conductivity, specific heat and 
density of fuel and cladding are adopted in the analysis. 

The manuscript is organized as follows. The governing equation 
of one-dimensional melting of a nuclear fuel rod and the corre- 
sponding enthalpy formulation of the problem are presented in the 
next section. The FD schemes adopted to the enthalpy formulation 
is described in the following section. The melting process of the fuel 
rod during full power operation are then presented. Subsequently, 
the effect of heat transfer coefficient between coolant and fuel rod 
on the temperature histories of the fuel pellet and cladding is 
investigated. The melting process in the nuclear rod with decay 
heat power after shutdown is also simulated. Finally, concluding 
remarks and some perspectives are provided. 

2. The mathematical formulation 

2.1. Transient heat conduction 

Consider the transient heat conduction in a nuclear fuel rod, 
which consists of heat-generating cylindrical pellets of uranium 
dioxide (UO 2 ) fuel sealed within a Zircaloy cladding tube, as shown 
in Fig. 1. A gap between the pellets and the cladding tube is filled 
with pure helium as a high thermal conductivity fill gas. We 
simplify the problem by assuming axisymmetry of temperature 
distribution and neglecting the axial heat conduction term and the 
spatial variation of the heat generation across the fuel rod. The 
temperature-dependent thermophysical properties are considered 
in the fuel and the cladding. Under the above assumptions, we have 
the following governing equations with appropriate boundary and 
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Fig. 1. Sketch of a nuclear fuel rod. 

initial conditions for one-dimensional transient heat conduction in 
a fuel rod [4,5]: 


The heat conduction equations (based on Eq. (1)) governing the 
temperatures in liquid and solid phases of fuel are given by the 
following: 

T ¥, hMw) +*• »<’■<*.«■ 

(7) 

Pfs (TVs) c Pk (r 6 ) ^ \ 1 (rk fs (r fs ) +g, s f (t) < r < R { , 

( 8 ) 


with the boundary conditions (based on Eqs. (2) and (3)): 



r = 0, 


( 9 ) 



dTk 

dr 





r = Rr 


(10) 


where the subscripts 1 and s refer to liquid and solid phases, 
respectively, Sf is the location of the solid-liquid interface in fuel 
which is not known a priori. 

The heat conduction equations (based on Eq. (4)) for cladding 
are given by: 

Pci(Td)Cp C i(Tci) = \ | (rkdOd) , R d< r < s c (t), 

HD 


o,(r,)c K (T f )f 


§ - 0 ’ r -° ■ 



0 < r < R fl 


(1) 

( 2 ) 


Pcs(Tcs)Cp cs (Tcs) = 7 I (rkcsiTcs) ) , s c (t) <r<R C0 , 

( 12 ) 

with the boundary conditions (based on Eqs. (5) and (6)): 


M t ’)w - R i h ‘( T '- T LJ- r - R - (3) 

Mr)c p .(Tc)^f = l§:(rW^, Rc,<r<Rc., (4) 

= r = R «- < 5 > 

r = R c „ (6) 


where 7> and T c are temperatures in fuel and cladding, respectively; 
Pf and p s , their densities; c Pf ,c Pc , the specific heats; kf, k c , the 
respective thermal conductivities; g, the volumetric heat genera- 
tion rate in fuel; fr g , the convective heat transfer coefficient for the 
gap, hoc the convective heat transfer coefficient between the clad- 
ding and the coolant, while Rf, R Cil R Co are the radius of fuel, the 
inner radius and the outer radius of cladding, respectively. R g is the 
mean of R Ci and R Co ; Too is the temperature of coolant. 

2.2. Transient heat conduction with phase change 


U-Tc), r-R ci , (13) 

-kcs(T^ = boo (Tcs -Too), r = R co , (14) 

where s c is the location of the solid-liquid interface in cladding. 

By performing an energy balance between the solid and liquid 
phase along the phase change front, we have the following inter- 
face heat flux equations: 

kn(Ta) d ^ + P fs (T fs y { ^ = fc fs (r fs )^, at r = s f (t), 

(15) 

Mr c i)^ + Pc S (Tc S )Lc^ = MTcs)^, at r = s c (t), 

(16) 

for fuel and cladding, respectively, where U and L c are the latent 
heats of melting. 

The temperature at the interface is the constant melting tem- 
perature, giving 


During a loss-of-coolant accident (LOCA), the initial stored en- 
ergy from operation and heat generated by the radioactive decay of 
fission products cannot be removed by the coolant. As a result of 
severe overheating, nuclear fuel rod meltdown occurs when a nu- 
clear fuel element exceeds its melting point. In this case, the region 
of fuel or cladding is divided into two regions: the liquid phase and 
the solid phase. 




— ^fim 


(17) 


7cllr=s c — ^cs| r =s c — ^ cm ’ 0^) 

for fuel and cladding, respectively, where Tf m and T cm are the 
melting temperatures. 
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2.3. Enthalpy formulation 


In Section 2.2, the energy equations are written separately for 
the solid and liquid phases, and the temperatures are coupled 
through the interface energy balance condition, which brings a 
large proportion of the available numerical schemes (e.g. finite- 
difference method) difficult to implement. To overcome 
this problem, the enthalpy form of the energy equations is 
adopted to generate a single equation applicable in both the 
phases. For fuel, Eqs. (7) and (8) can be written in the enthalpy 
formulation as 



1 

r dr 


(rkf(l>) 





0 < r < Rf, 


(19) 


rTf 

where H f = / c Pf (T f )dT f + 0L f according to Ayasoufi et al. [23] 
Jo 

and 0 equals 1 for the liquid phase and 0 for the solid phase. 
Assuming a constant value for the heat capacity inside each phase, 
we can obtain 


c Pfs Lf 7 Lf — Lf m , 

Cpfs^fm + Lf + c Pn (j'f ~ Tfm) 5 for Lf > Lf m . 


( 20 ) 


Equivalently, Eq. (20) can be written as 

(Hf/c Pfs , for H { <c p J fm , 

Lf == \ Tf m , for c Pfs Tf m < Hf < c Pfs Tf m + Lf , 

{ (Hf - c Pfs T fm - Lf) j c p „ + T fm , for H f > c Pfs T fm + L f . 

(21) 

For cladding, Eqs. (11) and (12) can be written in the enthalpy 
formulation as 


0H C i a / . _ . ar c 
Pc lT - 7 8? * (Iy ^F 


-g, 


R C i <r < R C{ 


( 22 ) 


where H c = [ c Pc (T c )dT c + 0L C . Assuming a constant value for 
the heat capacity inside each phase, we can obtain 


Hfi +1 = H fi + ( T f(i+1) - T fi) 

- ^ f^ 2 ) (b" - Vd) + G ’ '' • 12 N f-r 

(25) 

where r,- = iAr, <r+ = k f(1+1/2) At/(Ar 2 /) fi ), ^ = k f{i _ 1/2) At/(Ar 2 /) fi ), 
G=gktlpfi. The integers z and n increase monotonically along the r 
and t coordinates, respectively. Nf is the node index at the boundary 
of fuel. Since the thermal conductivities may be different at the 
positions z + 1/2 and i - 1/2 due to the phase change at node z, they 
should be ascertained for each node at each time step as follows: 


! ^f (i-l/2) — kf(i+l/2) — kfs? for Hfj <H f s ( max ), 

^f(i— 1/2) = kfl an d fcf(i+l/2) = fcfs* for ^fs(max) < H fi < H fl ( m in ), 
^f(i-l/2) — ^f(i+l/2) ~ ^fb for Hfi 

(26) 

where Hf S ( max ) and Hfi( m j n ) represent the maximum enthalpy value 
of solid phase and the minimum enthalpy value of liquid phase of 
the fuel, respectively. 

Similarly, the finite-difference form for the heat conduction 
equation of cladding (Eq. (22)) can be written as 

HcT 1 = H d + (“1^) ( T c(i+V ~ b"i) 

- ff ci (^=P) (IS - T cVd) > f = N ci + 1 A + 2 , . . „JVco - 1 , 

(27) 


where (T+j <= k C (i+i/ 2 )At/(Ar 2 p d ), ^ = k c{i _ 1/2) At/(Ar 2 /) ci ). N ci and 
N co are node indices at the inner and outer boundaries of cladding, 
respectively. When considering phase change, the thermal con- 
ductivities for each node at each time step can be expressed as: 

! ^c(i- 1/2) — ^c(i+l/2) = ^CS} for H c [ < H cs ( max ) , 

^c(i-l/2) — K\ a nd/< c(l+ i /2) = l<c s, for ^cs(max) — Lf fi <^cl(min)? 
^c(i-l/2) — ^c(i+l/2) — ^cb for W c j > , 

(28) 


( c p cs Tc, for T c < T C m, 

\ ^Pcs^cm + Lc + c pc i(T c — L C m), for T c > L C m- 


(23) 


Equivalently, Eq. (23) can be written as 

( H c /c Pcs , for H c <Cp cs T c m , 

T c = J ^cm, for C Pcs Tcm < H c < Cp cs T c m + L c , 

1 (^C - Cpcs^cm — Lc) j Cpcl + L cm , for He > C Pcs T C m + L c . 

(24) 


3. Solution method 

3.1. Explicit finite-difference scheme 

The enthalpy formulations of the melting problem of fuel rod 
are discretized by the explicit FD scheme. The “r— t" domain is 
subdivided into small intervals of constant Ar in space and variable 
At in time. For the heat conduction equation of fuel (Eq. (19)), the 
forward and central differential schemes are adopted for time and 
space derivatives, respectively, yielding the following finite- 
difference equation: 


where H CS ( max ) and H c i(min) represent the maximum enthalpy value 
of solid phase and the minimum enthalpy value of liquid phase of 
the cladding, respectively. 

Now we consider the finite-difference formulations for the 
boundary nodes. For the central point of the fuel (z = 0), we per- 
formed the integration of Eq. (19) as 


Ar/2 


Ar/2 


Ar/2 


/ A 9 ? 2lrrdr = / 7 J (rk f ^pj2nrdr+ J g2nr 


dr 


(29) 


By using the boundary condition, Eq. (2), the above equation 
gives 


Ar dHf dTf Ar 


6r 4 


(30) 


Then, discretize the equation with the forward scheme in time 
and the central difference in space, yielding: 

Hfo +1 = Hfo + 44 (b " 1 - bo) + G, (31) 
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where = kf(i/ 2 )At/(Ar 2 p f0 ), G = gAt/p/p. The thermal conduc- 
tivity /<f(i/ 2 ) can be calculated by Eq. (26). For the other boundary 
nodes, the difference formulation Eq. (25) is adopted for i = Nf, 
while Eq. (27) is employed for i = N C i and i = N co . In this circum- 
stance, some fictitious nodes (i = Nf + 1, N c i - 1 and N co + 1) are 
introduced into the formulations, which can be eliminated by dis- 
cretizing the boundary conditions Eqs. (3), (5) and (6), respectively: 


r rn 

f(JVf+l) 

1 

l 

II 

2A rRghg 
kfNfRf 

te 

~ T c\ ,), 

(32) 

^c(N ci — 1) 

— T n 4_ 

- J c(N ci +l) + 

2ArR g hg 
kcN a Rci 

( T m 


(33) 



2Arhoo / 

Ya, 

- T ) 

(34) 

J c(N C0 + 1) 

— 1 c(N C0 -l) - 

kcN co ' 

/cJVco 

J 00 J • 


<t>i = (Hfm - H£-)/(h£ +1 - H 'fi) for fuel melting, 

<t>i = (H cm - H( ! i )/(H( , + 1 - H” j for cladding melting. 

(40) 

At time t 2 - the temperature at node z is Tf m (or T cm ). The tem- 
perature at the other node points are easily estimated by using a 
linear interpolation: 

C = x ‘-( T fP +1 - T f P ) + r fP’ p for fuel ' 

T'cp X ‘ = Xi (r^ - T" p J + T" p , p*i, , for cladding. 

(41) 

4. Results and discussions 


The thermal conductivity k fN[ can be calculated as: 

J l<fN, = '<fs for H fN[ < Hf m , 

\ = k n for H m > H fm , 


where Hf m = Cpf S Tf m + Lf/2. Similarly, the thermal conductivities 
k cNc . and k cNco can be calculated as 


f ^cN ci — ^cs for H c n c . < H cm , 
1 l<cN a = fcd for H cNa > Hcm, 


and 


f kcN co — kes for H c n C0 < Hcm, 
\ k C N co — k c i for H cNco > H cm , 


respectively, where H cm = c pcs T cm + L c /2. 

Once the enthalpy value of each node at a given time step is 
available, the corresponding temperatures of fuel and cladding are 
determined from Eqs. (21) and (24), respectively. The nuclear fuel 
rod is initially in the solid state and at a uniform temperature: 

'T° = T 0 , for i = 0, 1, 

T° = 7 o, for i = N ci ,N ci + 1, ...,N C0 , nfn 

H° mc pfs T°, for imOJ,...,N f , 1 j 

. H° a = Cp CS T pj , for i = N ch N d + 1, „„N C0 . 

Note that the conditions <r£ , , a't, < 0.5 should be satis- 
fied to guarantee the stability of the explicit finite-difference 
scheme. 


3.2. Moving boundary 

When the phase change boundary is on the node i, the nodal 
enthalpy is the sum of the sensible heat of the phase change and 
half the latent heat associated with the phase change [ 1JJ, viz. 
Hfi = Hf m = Cpf s Tf m + Lf/2 or H c i = H cm = Cp CS T cm + L c /2. Therefore, 
whenever the enthalpy at a node point is such that at time level 
n + 1, HQ < Cpf S Tf m + Lf/2 (or H^- < Cp C sT C m +L c /2) and 

Hfi +1 > Cp fs T fm + Lf/2 (or H"-’ > c pC sT C m + L c /2), the phase change 
boundary has passed through the point zAx in the preceding time 
interval At. Furthermore, assuming that the enthalpy changes lin- 
early in any time interval, the time at which the phase change is on 
the point zAx is given by 


ti = (n + 0,)At, (39) 

where </>* < 1 is estimated via linear interpolation in time, viz. 


The geometrical parameters and thermophysical properties of a 
typical fuel rod in a PWR are presented in Table 1 [24]. The prop- 
erties of Zircaloy are used in place of Zircaloy-4 employed in the 
following calculation as the properties Zircaloy-4 are not available 
to us. According to IAEA [25] and Fink [261 , the material properties 
of fuel and cladding are presented in Appendices A and B, respec- 
tively, while the average over temperature of those properties for 
both solid and liquid phases are listed in Table 2. After performing 
the mesh convergence study, 20 and 5 nodes are, respectively, used 
for discretizing the fuel and cladding in the FD scheme. 

4.1. Melting under full reactor power 

Before analyzing the temperature profile during the LOCA ac- 
cident, a normal operating condition is simulated with hoc = 34 kW / 
m 2 K, considering as the initial condition before the occurrence of 
accident. The temperature history of the nuclear fuel rod at the 
positions of r = 0, ftf, R C [ and R co during the heat-up operation is 
shown in Fig. 2, where it can be seen that the steady-state tem- 
peratures achieve at around t = 40 s. Figs. 3 and 4 show the 
enthalpy and temperature histories of the nuclear fuel rod with 
constant material thermophysical properties under an accident 
condition (the heat transfer coefficient between coolant and fuel 
rod h oo =0.2 kW/m 2 K) associated with full reactor power opera- 
tion, respectively. From Fig. 3, jump discontinuities can be clearly 
seen, which prove the occurrence of phase change of cladding (r = 0 
and ftf) and fuel (R C i and R co ). In addition, we can observe that the 
variation of enthalpy and temperature across the cladding is small 


Table 1 

Parameters of a PWR reactor [24]. 



h g (kW/m 2 K) 


6.4695 



q' (kW/m 3 ) 


318,121 



Rf (mm) 


4.1 



R ci (mm) 


4.18 



Rco (mm) 


4.75 



Too (K) 


578 


Table 2 





Average properties of the fuel rod. 




Properties 

Fuel 


Cladding 



Solid 

Liquid 

Solid 

Liquid 

P (kg/m 3 ) 

10,330 

8220 

6415 

5878 

k (W/m K) 

2.8 

2.5 

30.8 

37.0 

c P (J/kg I<) 

421 

351 

390 

546 

tm (kj/kg) 

260 


155 


T m (K) 

3120 


2150 
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Fig. 2. Temperature history of the nuclear fuel rod during a heat-up operation 
(hoo = 34 kW/m 2 K). 


due to the higher thermal conductivity and shorter radial length 
compared to those of the fuel. The enthalpy method is also carried 
out with variable temperature-dependent properties to obtain the 
temperature history, as shown in Fig. 5. The solution requires a 
higher computational cost than that with constant average prop- 
erties. Figs. 4 and 5 show the temperatures at r = 0, Rf, R c \ and R co 
nearly reach the steady-state at t = 200 s. 

The whole analysis can be subdivided into the following four 
different phases: (1) transient phase before the clad melting, (2) 
clad melting phase, (3) transient phase before the fuel melting and 
(4) fuel melting phase. Table 3 gives the durations for the four 



Fig. 3. Enthalpy history of the nuclear fuel rod with constant material properties under 
an accident condition (h „ = 0.2 kW/m 2 K) associated with full power operation. 



Fig. 4. Temperature history of the nuclear fuel rod with constant material properties 
under an accident condition (h^ = 0.2 kW/m 2 K) associated with full power operation. 

phases of the nuclear fuel rod with constant and variable properties 
under the accident condition associated with full power operation. 
It can be observed that the relative error reaches up to 25.34% at the 
second phase and drops to 12.78% at the fourth phase. The higher 
errors for the first and second phases are due to the large variation 
of thermal conductivity for the solid phases of both cladding and 
fuel. Since the properties have a minor variation to the liquid 
phases, the lower errors are obtained for the third and fourth 
phases. In general, the solution with constant properties is over- 
estimated. Fig. 6 shows the temperature profiles along the radial 
direction of the nuclear fuel rod with constant properties when 
cladding completely melts (the moment at the end of the second 



Fig. 5. Temperature history of the nuclear fuel rod with variable material properties 
under an accident condition (h^ = 0.2 kW/m 2 K) associated with full power operation. 
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Table 3 

Durations for the four different phases of the nuclear fuel rod with constant and 
variable material properties under an accident condition (/!«,= 0.2 l<W/m 2 K) 
associated with full power operation. 


Properties 

Duration (s) 



Phase 1 

Phase 2 

Phase 3 

Phase 4 

Constant 

33.60 

4.65 

24.47 

67.94 

Variable 

29.52 

3.71 

27.96 

60.24 

Relative error 

13.82% 

25.34% 

12.48% 

12.78% 



(a) 



(b) 


phase) and fuel completely melts (the moment at the end of the 
fourth phase), respectively. Due to the small thickness of cladding, 
the temperature difference for fuel at r = 0 and r = Rf is consider- 
ably greater than that for cladding at r = R C [ and r = R co . 

Fig. 7 presents the transient phase change interface position of 
cladding and fuel with constant properties under an accident 
condition. For fuel, the interface initially moves relatively fast, then 



(a) 



(b) 


Fig. 6. Temperature profiles along the radial direction of the nuclear fuel rod when (a) 
cladding completely melts, and (b) fuel completely melts. 


Fig. 7. Transient phase change interface position of (a) cladding and (b) fuel under an 
accident condition (/!<* = 0.2 kW/m 2 K) associated with full power operation. 
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the velocity becomes lower, which is due to the existence of heat 
generation only in the region of fuel. 

4.2. Effects of heat transfer coefficient h ^ 

In this case, it is assumed that the reactor is not shutdown 
properly and is still generating full power when the coolant is lost. 
To investigate the effect of the heat transfer coefficient between the 
coolant and fuel rod, h the temperature histories of the nuclear 
fuel rod with constant properties at the positions of r = 0, ftf, R c \ and 
R co are calculated with hoc = 0, 0.01, 0.1 and 0.2 kW/m 2 K, as shown 
in Fig. 8. It can be clearly seen that the velocity of temperature 
increase of the fuel rod becomes gradually high with decreased 
heat transfer coefficient h 

The durations for each phase of the nuclear fuel rod under 
different accident conditions associated with full reactor power 
generation are presented in Table 4, where the solutions using the 
classical lumped model proposed by Chen et al. [2] are also given. The 
situation with hoo = 0 represents the worst loss-of-coolant accident, 
for which the duration of each step is shorter than that for other 



(a) 



(c) 


situations with h = 0.01, 0.1 and 0.2 due to the insulation boundary 
condition. The solutions based on the enthalpy method have good 
agreement with the ones provided by Chen et al. £2J, especially for 
the first step, hoo = 0.01 does not provide a significant change in time 
history of r = 0, Rf, R c \ and R co . As hoo increases, the durations of each 
step extend, while the classical lumped model fails to predict the 
exact melting behaviour of nuclear rod except the first step. 

4.3. Melting under decay power after reactor shutdown 

In order to simulate a possible real situation, such as the accident 
of the Fukushima nuclear power plant, we consider the situation in 
which the reactor is shutdown but is not cooled adequately. The 
decay heat generation rate can be expressed as a function of time 
after reactor shutdown and reactor operating time before shutdown 
[241 . In this case, we hypothesize that the accident happens one year 
after the startup of the reactor. Since the decay power initially de- 
creases to around 6% of the full power in normal operation, the 
melting of the fuel rod does not occur due to lower heat generation 
rate in the fuel rod and high heat transfer coefficient (such as 



(b) 



(d) 


Fig. 8. Temperature history of the nuclear fuel rod at (a) r = 0, (b) r = R f , (c) r = R ci and (d) r = R co under different accident conditions ( h „ =0, 0.01, 0.1 and 0.2 kW/m 2 K). 
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Table 4 

Durations for each phase of the nuclear fuel rod under different accident conditions 
(ha, =0, 0.01, 0.1 and 0.2 kW/m 2 K) associated with full power operation. 


Properties 

Duration (s) 




Phase 1 

Phase 2 

Phase 3 

Phase 4 

/loo =0 

Enthalpy method 

20.15 

2.47 

13.39 

9.97 

Classical lumped model [2] 

20.12 

1.85 

15.33 

8.49 

/loo = 0.01 

Enthalpy method 

20.57 

2.53 

13.61 

11.04 

Classical lumped model [2] 

20.47 

1.89 

15.86 

8.87 

/loo =0.1 

Enthalpy method 

25.27 

3.19 

16.47 

21.92 

Classical lumped model [2] 

/loo =0.2 

24.47 

2.30 

23.37 

14.93 

Enthalpy method 

33.60 

4.65 

24.47 

67.94 

Classical lumped model [2] 

31.64 

3.17 

53.85 

56.29 


h oo = 0.1 or 0.2). To force the melting process to occur, the core is 
considered as fully insulated, which means that there is no heat 
transfer between the fuel rod and the coolant. Fig. 9 displays the 
temperature history of the nuclear fuel rod with constant material 
properties under an adiabatic condition (hoc =0) associated with 
decay reactor power. It can be observed that at the initial stage, the 
temperature at fuel center (r = 0) drops sharply due to the reduced 
heat generation rate and the heat transfer coefficient. Besides, the 
temperature histories of each position are quite coincident due to 
the fully insulated condition. Table 5 shows the durations for each 
phase of the nuclear fuel rod associated with decay reactor power 
and the relative errors between the solution using the classical 
lumped model [2J and the solution using the enthalpy method 
proposed in this work. The greatest relative error is 6.58%, con- 
firming the applicability of the classical lumped model for this case. 

5. Conclusions 


In the present study, the enthalpy method with explicit finite- 
difference discretization is adopted to analyze the melting process 



Fig. 9. Temperature history of the nuclear fuel rod with constant material properties 
under an accident condition (h x = 0 kW/m 2 K) associated with decay heat power. 


Table 5 

Durations for each phase of the nuclear fuel rod under an accident condition 
(h oo = 0 kW/m 2 I<) associated with decay heat power. 


Properties 

Duration (s) 



Phase 1 

Phase 2 

Phase 3 

Phase 4 

Enthalpy method 

655.5 

44.1 

778.4 

447.6 

Classical lumped model [2] 

673.7 

41.2 

786.5 

447.5 

Relative error 

2.78% 

6.58% 

1.04% 

0.02% 


of a nuclear fuel rod under the accident conditions of abnormal 
cooling. The heat transfer coefficient between cladding and 
coolant and the volumetric heat generation are two predominant 
factors in determining the temperature histories and phase change 
interface positions of fuel and cladding. To reduce the computa- 
tional cost, the average material thermophysical properties are 
employed, and the relative error between the solution with con- 
stant properties and the one with variable properties reaches up to 
25.34% at the second phase and drops to 12.78% at the fourth 
phase. Due to the small thickness of cladding, the temperature 
difference for fuel at r = 0 and r = Rf is considerably greater than 
that for cladding at r = R C [ and f — Rco- Through parametric study, it 
was found that the velocity of temperature increase of the fuel rod 
becomes gradually high with decreased h oo. By considering the 
shutdown reactor with decay heat power, we observe that the 
temperature histories of r = 0, Rf, R c i and R co are quite coincident 
due to the fully insulated condition h x = 0. Besides, it can be also 
seen that the classical lumped model cannot give the accurate 
results when the heat transfer coefficient becomes higher, there- 
fore, the lumped model should be further developed to enhance its 
applicability. 
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Appendix A. Material properties of fuel (UO 2 ) 

The thermophysical properties of UO 2 are described by the In- 
ternational Atomic Energy Agency (IAEA) [25]. The expressions for 
enthalpy of fuel are (unit: kj/mol): 

H ft (Tf) = -21.1762 + 52. 1743t + 43.9753t 2 -28.0804I 3 
+ 7.88552t 4 - 0.52668t 5 
+ 0.71391t-\ for 298.15 I< 

< Tf < 3120 K (solid phase), 

(A.1) 

where t = Tf/1000, and 

Hfl(Tf) = 8.0383 x 10 2 +2.5136 x 10“ 4 T f - 1.3288 
xl0 6 /r f , for 3120 K 

<T f <4500 K (liquid phase). (A.2) 

Note that 1 mol of UO 2 is about 0.270 kg. 

The expressions for specific heat of fuel are as follows (unit: J/ 
mol K): 
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Appendix B. Material properties of cladding (Zircaloy) 

CpfjTf) = 52.1743 + 87.951 t - 84.241 It 2 + 31.542t 3 

' ' According to IAEA [25], the expressions for enthalpy of fuel are 

-2.6334t 4 -0.71391 t- 2 , for 298.15 K (unit: kj/mol): 

< T f < 3120 I< (solid phase), (A.3) 


Hcs(Tc) 


-7827.595 + 24.16187c + 4.37791 x 10- 3 T 2 + 6.9942 x lO 4 !” 1 , for 298.15 K< 7 C < 1100K (solid phase), 
-525.539 + 25.60747 c + 3.4008 x 10- 4 T 2 + 1.9458 x 10 3 7 3 + 2.2843 x 10- 10 7 4 + 5.0466 x 10 4 7- 1 

for 1100 I< < 7 C < 2150 K (solid phase), (B.l) 


c pf i(Tf) = 0.25136 + 1.3288 x 10 9 /7 f 2 , for 3120 K < 7 f 
<4500 K (liquid phase). 

(A.4) 

The expressions for thermal conductivity of fuel are (unit: W / 
m K): 


H d (7 c ) = 768.563 +52.389 x 10- 3 7 c - 0.956331 x 10- 2 T 2 
+ 2.53317 x 10 6 7 3 - 1.52605 
x 10- 10 7)?, for 2150 K 
< 7 f < 4100 1< (liquid phase). 

(B.2) 


‘W 


100 


6.548 + 23.533t 


6400 

+ ^T ex P 


-16.35 


(A.5) 


for 298.15 K < 7 f < 3120 I< (solid phase), 


Note that 1 mol of Zircaloy is about 0.091 kg. 

The expressions for specific heat of cladding are as follows (unit: 
J/mol I<): 


{ 255.66 + 0.10247 c , for 298.15 K < 7 C < 1100 K (solid phase), 

255.66 + 0.10247c + K> for 1 100 I< < 7 C < 1214 K (solid phase), 

597.1 - 0.40887c + 1-565 x (10“ 4 )7 2 + ? for 1214 K < 7 C < 1320 I< (solid phase), 
597.1 - 0.40887c + 1.565 x (lQ- 4 )7 2 , for 1320 I< < 7 C < 2150 I< (solid phase), 


and 

kfl^Tf) = 2.5, for 3120 K < 7 f < 4500 I< (liquid phase). 

(A.6) 

The expressions for density of fuel are (unit: kg/m 3 ): 

P/s(7/) ='^f 63 . ^ 298.15 K<7 f 

< 3120 K (solid phase), (A.7) 

where 

? = 9.9734 x lO” 1 + 9.802 x 10-®7 f - 2.705 x 10- 10 7 f 2 + 4.391 x 
10 13 7 3 for 298.15 I< < 7 f < 923 K and 

? = 9.9672 x lO” 1 + 1.179 x 10- 5 7 f - 2.429 x 10- 9 7 f 2 + 1.219 x 
1 0 12 7 3 for 923 I< < 7 f < 3120 K, and 

Pfl(jf) = 8860 “ 0.9285 (r f - 3120), for 3120 K < 7 f 

< 4500 K (liquid phase). 


where ? = 1058.4 x exp[(7 c - 1213.8 ) 2 /719.61], and 


CpJTc) = (52.389 - 1.912661 x 10“ 2 7 c + 7.5995 

x 10 _6 7 2 - 6.1042 x 10 _1 °7 3 ), for 2150 K 
< 7 f < 4100 K (liquid phase). 

(B.4) 

The expressions for thermal conductivity of cladding are (unit: 
W/m K): 


kcs(Tc) = 12.767 - 5.4348 x 10~ 4 T C + 8.9818 
x 10- 6 T 2 , for 298.15 K 

< 7 C < 2150 1< (solid phase), (B.5) 


kc (7 C ) = 36.5, for 2150 K < 7 f < 4100 K (liquid phase). 

(B.6) 


(A.8) 


The expressions for density of cladding are (unit: kg/m 3 ): 
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/ 6595.2 - 0.1477Tc, for 298.15 K < T c < 1100 I< (solid phase), 
\ 6690 - 0.18557c, for 1100 K < 7 C < 2150 I< (solid phase), 


(B.7) 


p d (7 c ) = 6844.51 - 0.6098987 c + 2.05008 x 10“ 1 * * 4 * * * T C 2 
-4.47829 x 10“ 8 * * * 7 2 +3.26469 
x 10“ 12 7^, for 2150 K 

< T f < 4100 K (liquid phase). (B.8) 
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